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ABSTRACT 


Seismic reflection data recorded close to the 
geophysical observatory at Edmonton, Alberta have been 
interpreted to produce a crustal model for this area. 
Velocity analysis was performed with the aid of a computer 
program, written specifically for this work, which does not 
require common depth point data and which will allow for 
Substantial dip on reflecting interfaces. The model shows 
the presence of 15 degree southeasterly dips within the 
crust, with the base of the crust essentially flat. Total 


crustal thickness of 35.5 kilometers is indicated. 
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CHAPTER 1 


iNTRODUCTION 


1.1 History of Crustal Studies 


The earth's crust or outer shell has always been of 
greater interest to seismologists than its small size would 
suggest. There are several reasons for this. Since this 
region completely encloses the rest of the earth it will 
affect any measurements made on deeper regions. Without a 
good knowledge cf crustal structure it would be difficult to 
obtain precise information about the interior of the earth. 
The second reascn for the interest in the crust is the 
effect it has on man. The crust is the source of the 
earth's ores and minerals and a knowledge of the 
relationship of crustal properties to location and formation 
of these reserves is valuable both economically and to our 
weli-being. Another important consideration is the 
earthquake prediction problen. If we hope to predict 
earthquakes then a detailed knowledge of the crust and upper 


mantle may be important. 


The first studies of crustal structure were made 
through observations of near earthquakes. The initial 
discovery occurred in 1910 when Mohorovicic identified the 
velocity discontinuity which now bears his name. This is 


referred to as the M discontinuity and is generally accepted 
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as the base of the crust. Tne71925e7Conrad  ereported) .a 
discontinuity above the M which does not seem to be as 
universal as the M discontinuity. fhe Conrad discontinuity 
is generally regarded as the transition between a granitic 
crust above and an intermediate layer of more basic 


composition below. 


Explosion seismology has proven to be the most powerful 
tool for studying crustal structure. This method was first 
applied in petroleum exploration in 1921 by Drs. J. Clarence 
Karcher, William P. Haseman, Irving Perrine and Mr. William 
C. Kite (Schriever, 1952). None of these people could have 
predicted at that time the success this method would have. 
It is today the most widely used and effective of ail 
petroleum prospecting methods (Alien, 1971). These 
techniques have been applied to studies of the deep crust 
and both refraction, and more recently reflection, methods 


have been employed. 


The first published report of deep crustal reflections 
was made by Junger (1951). During normal seismic work in 
Big Horn County, Montana the recording camera was allowed to 
run for 10 seconds to study a long-lasting, low-frequency 
surface disturbance set up by a 25 pound charge. The 
resulting record showed a reflection at 8.5 seconds. Other 
records were obtained in adjoining areas which showed 
reflections at times of 7.0 to 7.5 seconds. The recording 


of these reflections was found to be repeatable, and on the 
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basis of the energy they contained, Junger argued that they 
could not be muitiples and must be reflections from inside 
the basement. Other isolated instances of deep reflections 
were reported and a good summary of these reports and the 
arguments for and against the existence of such reflections 


can be found in Steinhart and Meyer (1961). 


A later review by James and Steinhart (1966) describes 
the results of reflection and refraction crustal studies 
from 1960 to 1965. An excellent recent review of seismic 
crustal studies in all parts of the world can be found in 
Chapter 3 of the Geophysical Monograph 13 (Hart, 1969). 
This work contains a collection of papers by different 
authors pertaining to studies in many countries. A brief 


Summary of this work follows. 


Seismic explosion studies in western Europe have been 
made by reflection and refraction surveys of petroleun 
companies, quarry explosions of mining companies and 
scientific explosions in the sea, Alpine lakes and in 


boreholes (Closs, 1969). 


In the Alpine region, refraction work has been done 
mainly by workers in France, Germany, Italy and Switzerland. 
They found a general thickening of the intermediate layer 
along the southern inner margin of the Alps. A great deal 
of work was done in central France, also. Shots in the 
North Sea and the Irish Sea have led to a crustal model with 


a very thin, or non-existent intermediate layer. 
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In northern Europe surveys have been carried out in 
Finland, Sweden, Norway and Denmark. Good review articles 
for each of these countries can be found in a_ publication 
edited by Vogel (1971). Velocities and thicknesses of 
granitic and basaltic layers as well as mantle velocities 
are given for Finland, Sweden and for three profiles in 
Norway. A refraction study in Denmark has suggested a very 


thin granitic layer in that area. 


in Germany a great deal of crustal structure work, much 
of which is based on reflection seismology, has been done. 
Most of the German work is statistical in nature, involving 
histograms which show the number of reflections in a certain 
area as a function of reflection time. A comprehensive 
review of the German work on deep crustal reflections is 
given by Dohr and Fuchs (1967). They found that deep 
reflections only correlate over short distances and exhibit 
considerable scatter in reflection times. A short review 
article which presents a diagram of the deep crustal 
structure of Germany and provides an extensive list of 


references of German work is given by Stein (1971). 


An analysis of wide-angle reflections in Germany and 
the USSR led Meisner (1967) to suggest that the M 
discontinuity is stepwise, interrupted by rhythmically 
arranged series of partial melts of lower velocity. Fuchs 
(1969) demonstrated that a simple layered model of the 


reflecting horizons in the earth's crust was not consistent 
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with certain observations, especially the large amplitudes 
and low cut-off frequency of deep reflections. He suggested 
a laminated transition zone with a series of velocity 


reversals as a model for deep crustal reflectors. 


Crustal studies in southeastern Europe have been 
carried out by several countries and are reviewed by 
Sollogub (1969). Most results come from refraction 
profiles, although good reflections from the M discontinuity 
were observed in Hungary. The earth's crust in the 
Carpatho-Balkan area was found to be a tlayered structure 
with seismic interfaces in addition to the M and Conrad 
discontinuities. The existence of mountain roots under the 
Crimean structures, Carpathians and Dinarides was proven and 
evidence of “anti-roots" along the Conrad discontinuity was 
found. MThicknesses and velocities of individual layers were 


mapped throughout these countries. 


In Russia results have been obtained with the deep 
seismic sounding (DSS) method designed by G. A. Gamburtsev 
in 1948 to 1955 and commenced in 1956 (Kosminkaya et al, 
1969). A description of the Soviet techniques and 
observational methods along with information about use of 
reflected and refracted waves can be found in Kosminkaya and 
Riznechenko (1964). This review also gives cross-sections 
and a discussion of the nature of stratification. Another 
excellent review which describes almost all aspects of the 


DSS program up to 1962 is given in a collection of Soviet 
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papers edited by Zverov (1967). 


Kosminkaya et al 


different tectonic zones within the USSR. 


(1969) present crustal models of 17 


These models give 


layer thicknesses and boundary velocities of the major 
eraustal layers. They are all characterized by block 
structure, different depths to the M discontinuity, 


different layer thicknesses, 


an M boundary velocity of 8.0 


The Vela Uniform program 
provided a 
America. 
and velocities 


layer were 


Many of these results are reviewed Healy and Warren 


The 


range of 5.9 to 6.2 km/sec and the 


ooG €toe8/. 15 knysec* In 


layer was found with a veiocity of 7.1 to 7.4 


P-wave velocity 


8.3 km/sec with the lower 


southwestern United States. 


Although much refraction 


PELOESHte 4149697 


Refraction studies to determine crustal 


very few reflection studies were made. 


different layer velocities and 


to 8.2 km/sec. 


and the Upper Mantle Project 


strong impetus to expiosion seismology in North 


thickness 
made on about fifty profiles. 


(1969). 


upper crustal layer was found to have a P-wave velocity 


lower layer a range of 


some areas a second deep crustal 


km/sec. The 


of the uppermost mantie ranges from 7.8 to 


velocities predominant Bt 


work was done in North America 


In 


fact, prior to 1968 the only North American group to use the 


reflection technique for studying crustal structure 


the University of Alberta (Clowes, 


was at 


1969). This work will be 


reviewed later as it forms a prelude to this thesis. 
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Studies of crustal structure of the Atlantic Ocean 
began shortly after World War II (Ewing, 1969). By 1960 
over one hundred refraction profiles had been recorded. One 
of the most significant results was the discovery that the 
crust is much thinner under the oceans than under the 
continents. The typical crustal structure consists of a 
sedimentary zone at the top with a P-velocity of 1.3 to 3 
km/sec, a second layer of velocity 4.5 to 5.5 km/sec anda 
lower layer of velocity 6.5 to 7.0 km/sec. Mantle P-wave 
velocity ranges from 7.4 to 8.4 km/sec beneath the ocean. 
The M discontinuity has an average depth under the Atlantic 
of 12 km and apparently disappears under the mid-Atlantic 
ridge. Similar structural features were found under the 


Pacific and Indian Oceans (Shor and Raitt, 1969). 


Japan too, has been active in seismological studies of 
the crust. A model of crustal structure for a cross-section 
of north-eastern Honshu has been prepared (Research Group 


for Explosion Seismology, 1968). 


1.2 Crustai Studies at the University of Alberta 


The PLrst seismological determination of crustal 
thickness in Alberta was that of Richards and Walker (1959). 
A team of fifteen seismic parties made measurements along an 
81 mile refraction profile. This profile was parallel to 
and about sixty miles east of the frontal thrust of the 


Rocky Mountains. They obtained a depth of 43 km to the 4 
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discontinuity and an upper mantle velocity of 8.2 km/sec. 
In addition they identified the Conrad discontinuity dipping 
at about 3 degrees to the south at a depth of 29 kilometers 
and with a velocity of 7.2 km/sec. Some evidence was found 


for other additional layers. 


The first use of explosion seismology at the University 
of Alberta for the study of crustal structure was made by 
Cumming, Garland and Vozoff (1962). Several shots were made 
at each end of an east-west refraction profile 250 
kilometers long in southern Alberta. They arrived at a 
preferred model which consisted of three layers between the 
top of the Precambrian and the Mantle. The uppermost layer 
was a Precambrian basement basic phase with a mean thickness 
of 20 km and a velocity of 6.4 km/sec. A velocity reversal 
was postulated and the thickness of a layer with velocity of 
6.1 km/sec was 11 km. A deep crustal layer with a veiocity 
of 7.32 km/sec overlay the mantle which had a mean depth of 
47.5 km and velocity of 8.25 km/sec. In all cases dip was 


less than one degree. 


Cumming and Kanasewich (1966) reported the results of 
an extension of the previous work along two reversed 
refraction profiles extending east and west from Suffield, 
Alberta for a total distance of 500 km. In addition, a one- 
way refraction profile extending west into the Rockies was 
shot. The revised crustal model which they arrived at from 


these studies consisted of three layers between the basement 
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and upper mantle. The uppermost layer referred to as the 
basement layer was found to have a velocity of 6.1 km/sec. 
Thickness of this layer ranged from 15 kilometers in some 
areas to zero just west of Suffield. The disappearance of 
this layer correlated with a steep gradient in the Bouguer 
gravity anomaly along the seismic line. The second layer, 
referred to as the sub-basement, extended to depths in the 
Eangereot 33 ttesg 37 km and hadea velocity of 6.5,ku/sec. 
Below this, the third or intermediate layer had a velocity 
of 7.15 to 7.2 km/sec and was underlain by the upper mantle 
with a velocity of 8.1 to 8.3 km/sec. The total crustal 


thickness ranged from 43 to 49 km. 


The first observation of deep reflections in Alberta 
was made by Robertson (1963) .. During the course of routine 
seismic reflection work in southwestern Alberta deep 
reflections at times of three to five seconds were observed. 
These were shown to be reflections from within the 
Precambrian basement and not multiples. The reflecting 


layers had a dip of eight degrees or more. 


The first deep reflection work at the University of 
Alberta was a study of near-vertical-incidence seismic 
reflections at Lomond, Aiberta by Kanasewich and Cumming 
(1965). They recorded good quality reflections at a time of 
11.4 seconds along an expanding spread. The use of 
horizontal geophones showed that the event was near-normal 


incidence. An x-squared, t-squared analysis gave an average 
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velocity from the top of the Mississippian to this reflector 
of 6.37 km/sec and a depth to the reflector of 34 km. This 
WaS in good agreement with the depth to the discontinuity 
observed in the refraction profile reported by Cumming et al 
(1962). However, since the velocities were different than 
for the Conrad in Europe and in order to avoid implying 
intercontinental correlation, this discontinuity has been 
referred to locally as the Riel or R discontinuity as 
Suggested by Professor D. H. Hall of the University of 


Manitoba. 


AS a result of the success of the preceding work, the 
deep crustal reflection program at the University of Alberta 
was expanded. The results of these studies were first 
reported by Clowes, Kanasewich and Cumming (1968). Seismic 
reflections over four profiles for a total of about 99 
kilometers in southern Alberta were studied. Power spectral 
calculations showed that the energy of the reflected 
wavelets was concentrated in the five to fifteen Hertz 
range. Along one profile the reflection from the R 
discontinuity was correlated over nearly 25 kilometers. 
Analysis of this data gave a revised average vertical 
velocity of 6.2 km/sec to a depth of 34 kn. Cther 
reflections were recorded from layers near the base of the 
crust including a layer corresponding in depth with the M 
discontinuity as determined by refraction studies. 
Continuous profiling along a 40 km line allowed construction 


of a seismic cross section which was characterized by 8 km 
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of relief over a horizontal distance of 25 km. This pointed 
out how the reflection method could be used to nmap 
complicated structure involving steep dips within the deep 


GRUSst 


These deep structures have been interpreted as a 
Precambrian rift valley lying below the flat sediments 
(Kanasewich, 1968). Seismic velocity data combined with 
measured Bouguer gravity anomalies led to a crustal model 
indicating that this rift valley was filled with lower 
density material. This rift was traced by gravity and 
magnetic trends for several hundred kilometers across 


Alberta and into British Columbia (Kanasewich et al, 1969). 


All of these previous studies had been concerned with 
obtaining crustai modeis of velocity-depth structure. In 
1970 Clowes and Kanasewich studied the attenuation 
properties of the crust andthe nature of the reflecting 
horizons (Ciowes and Kanasewich, 1970). Using a synthetic 
seismogram program which included attenuation as function of 
frequency and depth, they attempted to cbtain a model of the 
Q structure of the crust. They suggested that Q varies 
considerably above the basement with an average of about 
300, and that below the basement Q is approximately 1500 and 
probably increases with depth. Infaddition: to O;studies 
they used spectral analysis to examine the nature of the 
transition zones. They concluded that first-order 


discontinuities, linear gradients and step-like velocity 


=< 
2. ie o-Ps 
besaiog sivT 22 io so0sfeth fesacsizot » tevo teak 


si ot beeu 7) vod bodiae noritosize,s pal — 


rs j yeeatve palrys { istoutse 5b 
i 
d Vi injovi3ja _assp seedt 
ried iisv ~-2 723 I 
P . 
¥ ae | a ‘ * 
; : ef. FHOCUSB VILVSEDD -: ee 1 GSitee 


SO  epeaevs as S4iv Juemeasd sAd Bvod ie 


bck ¢ “daympely oe volun read. 


a a 


pfltgeh- Ate 


i ip 


increases gave generally unfavourable spectral 
characteristics. They found that a transition zone model 
consisting of sills of alternating high and low velocity 
material was in agreement with the observed data. The 
layers in this model were less than .2 kilometers thick and 


extended over a zone less than 1 kilometer in thickness. 


Chandra and Cumming (1972) continued the interpretation 
of refraction data in the area. They reported the 
interpretation of a 760 km refraction profile from Greenbush 
Lake, British Columbia to Swift Current, Saskatchewan with 
all the clder refraction work included. They examined this 
seismic refraction data in conjunction with reflection 
results and gravity and magnetic surveys in order to arrive 
at a unified interpretation for the area. Upper crustal 
velocities «vary laterally from 6.1 to 6.5 km/sec and 
boundaries between the velocity zones are most probably near 
vertical faults involving the whole crustal section through 
to the M discontinuity. The areas of 6.5 km/sec upper 
crustal velocity are in regions where the 6.1 km/sec layer 
is absent. The velocity beneath the two layers is 7.15 to 
7.17 km/sec and the upper mantle velocity varies from 8.0 
km/sec west of the Rockies to 8.3 km/sec in a region between 
Vulcan and Suffield where the profile crosses the Sweetgrass 


Arch. 


In addition to the explosion seismic studies of the 


crust in Alberta some attempt has been made to determine 
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crustal structure by studying the P-coda on earthquake 
seismograms. Ellis and Basham (1968) compared vertical to 
radial spectral ratios with theoretical ratios calculated 
using the Thomson-Haskell matrix formulation. They noted 
large deviations between experimental and theoretical 
curves, particularly at the University of Alberta Edmonton 
seismological observatory, and they attributed these 
deviations to shear wave conversion and scattering in the 
crust. They also found that for 20 events the cbserved 
azimuth was approximately 18 degrees more northerly than the 
true azimuth. This could be explained by localized dips of 


approximately 15 degrees on the crustal boundaries. 


In a more recent paper Somerville and Ellis (1972) 
Suggest a crustal model for the same area which gives better 
agreement with experimental P-coda studies. This model was 
obtained by studying vertical to radiai spectral ratios with 
the aid of synthetic seismograms. They found that a 
reasonable comparison between theoretical and experimental 
results could be obtained by including a low velocity layer 


several kiiometers thick at a depth of 12 kilometers. 


Another application of the P-coda spectral ratio method 
to the problem of determining crustal structure was made in 
the same area by Sprenke (1972). He could find no exact 
match between a theoretical and an observed spectral ratio 
but he did suggest the pone Te aiey that the upper crust is 


typified by 1) a very low S-wave velocity in surficial 
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layers, 2) a Poisson's ratio of one-third in the softer 
sedimentary layers and 3) a low velocity layer approximately 
3.5 km thick in the basement immediately beneath the 


sediments. 
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CHAPTER 2 


THE PROJECT 


a. teLocroduction 


As a result of the success with deep crustal 
reflections in southern Alberta, it was decided that work 
Should be done near the seismic observatory at Edmonton. 
The interpretation of P-coda studies would be facilitated by 
a detailed model of the crust in this area. An additional 
consideration was the availability of well velocity surveys. 
Figure 2.1 shows the location of the site chosen and nearby 


wells for which velocity information was available. 


Figure 2.2 shows the actual layout of the profiles for 
the shots. The recording set-up consisted of twelve groups 
of geophones spaced with 880 feet between centres. This, 
along with the 440 feet of spacing at each end, gave a total 
layout of two miles. Each receiver group consisted of 
sixteen geophones LaLGesOut Fin) as tapered array, the 


advantages of which are described by Clowes (1969). 


Three spreads were used for the collection of data in 
the summer of 1971. On August 25 the recording array was 
along the profile AB in figure 2.2 and a 40# and two 20# 
shots were fired at point A. The profile was along CD on 
September 2 and 40# shots were fired one, two and three 


Miles from each end at points A, F2, P3, P4, P5 and P6é. 
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Figure 2.1 Project area map. 
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Figure 2.2 
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For the shots at A, P2, P4 and P5 a pattern of holes as 
discussed by Clowes (1969) was used. For the remote shots 
at P3 and P6é patterns were not deemed necessary. The third 
spread, EF, was used on September 23 for five shots in 
Single holes. One 5# shot was at A and four 10# shots 
Started at F with successive 440 foot displacements towards 


A. 


2.2 The Digital Recording System 


Figure 2.3 is a schematic of the circuit for each 
channel. The seismometer shown represents the output from 
the 16 geophone array of Electro-Tech EVS-4 7 1/2 Hertz 
seismic detectors. Each data channel consists of the output 
of one of these arrays. RL is the line resistance and is 
proportional to the distance from the array to the recording 
position which is in the center of the 12 channel spread. 
RX iS a variable resistance such that RL + RX = 700 ohms. 
This value of resistance provides the optimum damping of the 
geophones in crder to give flattest frequency response. 
Thus, the output of each channel is proportional to the 
factor R'*/(R'tRL) where 1/R*=1/RX+1/10k. Since RL varies 
with position of the array then each channel must be 
normalized by the appropriate factor (1+RL/RX+tRL/10k) before 


amplitudes can be compared. 


In addition to the twelve data channels two other 


channels were recorded, one for the shot instant from a 
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Schematic of seismometer line resistance. 
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radio receiver and the other for absolute time from WWVB. 
These fourteen channels were recorded digitally on a 9-track 
Synchronous tape transport with a tape velocity of 6.25 
inches per second. The recording density is 800 bpi and 


hence the data transfer rate is 5000 bytes per second. 


The signal is converted from analog to digital form by 
a converter, with a resolution of 13 bits plus sign, 
providing a dynamic range of 84 db. Two 8-bit bytes are 
required to store each sample of the input signal and hence 
the sampling rate is 2500 times per second. Since there are 
14 channels, each channel is sampled once every 14/2500 
(.0056) seconds. This results in a Nyquist frequency of 
89.29 Hertz. One of the two unused bits in each two byte 
sample is used tc mark the end of one data conversion frame 


for all channels. 


In writing the records onto tape a length of 8192 bytes 
(4096 samples) per channel was chosen. This yieids a time 
interval of 22.94 seconds. Since fourteen channels were 
being recorded the block length on tape is 114,688 bytes. 
The inter-record gap resuits in the loss of .096 seconds of 


data when more than one block is written. 


Figure 2.4 is a diagram of the complete system. fhe 
tape recorder is a Model 7830-9 write-only transport made by 
Peripheral Equipment Corporation. The multiplexer, analog 
to digital converter and 5 kHz crystal clock plus associated 


interfaces and ccentrols are contained in a modified Model 
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Figure 2.4 
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#120 Data Acquisition System which was manufactured by Datum 
Tne. Controls allow the selection of from one to twenty 
channels for recording purposes and allow the choice of 


record length up to the maximum of 4096 samples per channel. 


Figure 2.5 shows the amplifier circuit. An input 
transformer of 18 db gain drives a 40 db fixed gain pre- 
amplifier which feeds the main amplifier through a high-pass 
two pole Bessel filter. A Bessei filter was used because of 
its linear phase shift properties. The latter may be by- 
passed. The main amplifier has switchable gain and its 
output passes through a four pole Bessel aliasing filter 
with a 50 Hz corner. Figure 2.6 shows the amplifier gain 


when the high-pass filter is by-passed. 


Figure 2.7 illustrates the type of data obtainable with 
this systen. The upper halif of the figure is part of the 
photographic monitor record of a 40 1b shot recorded at 88 
db gain setting. The lower traces are a Cal-Comp plot of 


the digitized data. 


This system has been described by Allsopp et al (1972). 


2.3 Velocity Log above the Basement 


In order to oktain reasonable values for the velocity 
in the lower crust it is essential to have reliable 
information about the velocity distribution in the overlying 


rocks. Therefore continuous velocity logs were obtained for 
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Figure 2.6 
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digital 


records. 
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the Big Hay Lake and the Canberra Millet wells shown in 
figure 2.1. The original logs gave average travel times for 
ten foot intervals, but no total travel times to formations. 
These were simplified by combining all regions where 
interval velocities in successive 10 foot layers were within 
3 per cent of each other into one thicker layer, with travel 


time being conserved. 


Depths to formation tops were also available for all 
wells shown in figure 2.1. These showed the sedimentary 
layers in this area to be essentially flat. Close agreement 
exists in depths for the Big Hay Lake and lLooma wells. 
Since the Looma well was near to the reflection profile, it 
was decided to adjust the velocity log for Big Hay Lake to 
fit this well. Both times and depths were available for the 
Looma well. This adjusted log would be used as a velocity 


model for the sedimentary rocks in the area of the profile. 


In order to match the formation depths in the two wells 
it was necessary to correct for the 150 feet greater depth 
to the top of the limestone in the Big Hay Lake well. 
Electric logs were obtained for the two wells and, by 
comparison of these logs, it was possible to ascertain which 
sections of the upper sediments should be removed from the 
velocity log in order to make the Big Hay Lake velocity 
survey compatible with the Looma well. In addition to this 
it was necessary to change the Wabamun velocity from 18,300 


to 19,200 ft/sec in order to match travel times in the 
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Paleozoic section. 


Since the velocity log began at a depth of 570 feet, a 
velocity model based on the hammer seismic results reported 
for this area by Sprenke (1972) was chosen for the top of 
the log. Also, the continuous velocity log did not extend 
to the basement, and so it was necessary to construct a 
velocity model from the travel times given for the botton 
section of the Leduc #530 well which was drilled to the 
basement. The location of this well is shown in figure 2.1. 
Table 2.1 compares elevations relative to sea level and one- 
way travel times for the Looma well and the adjusted 
velocity log, and shows the good agreement which was 


obtained. 


A sample seismic record from near this area was 
obtained from Imperial Oil. This was compared to a 
Synthetic seismogram calculated using the adjusted velocity 
log. The comparison was quite good although we did not know 
the effect of automatic gain control or the nature of the 
filtering done on the Imperial record. It is felt that this 
procedure yields a good approximation to the velocity 
distribution in the sediments and will provide an adequate 
model on which to base further analysis of deep reflection 


data. 
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TABLE 2.1 

FORMATION ELEVATIONS (FT) TRAVEL TIME (SEC) 

Looma New Log Looma New Log 
Lea Park 670 670 - 200 - 200 
Second White Specks -480 -480 wo2e eel 
Viking 200 -900 ~266 - 364 
Wabamun eT he ate O 444 ~ 443 
Cooking Lake -3060 sth) PAS, wo2D 
Beaver Hill Lake =335'10 -3310 soot - 540 


pikT Point -4010 -4010 2576 2574 


Pe 


CHAPTER 3 


COMPUTER PROGRAMS 


3.1 Methods of Velocity Analysis 


Determining seismic velocities by means of reflection 
profiles is a very old technique in the fieid of geophysics. 
One of the first reports of the method was given by Green 
(1938). He described the classic method of plotting 
recording distance squared against reflection time squared 
to give a straight line with slope equal to velocity 
squared. By using this method it is possible to obtain 
velocities within 3 per cent of the correct values (Gardner, 
1947). A detailed discussion of this method was given by 
Dix (1955) < All of these papers referred to common depth 
point shots. More recently, particularly during the last 
three or four years, the seismic exploration industry has 
given considerable attention to the development of computer 
velocity analysis systems (Montalbketti, 1971). A brief 
discussion of some of the principles involved in these 


methods will be given. 


In all of these methods it is assumed that the common 
depth point method has been used in recording the data. 
This method was outlined by Mayne (1962). Additional 
discussions of this method were given by Mayne (1967) and 


Courtier and Mendenhall (1967). The method provides for 


. » ] a ar i 
a i] 7 9 , 7 
in 7 a 
7 = a 
- 
7 
S ASC? ae 
- ine: 4HO 
J 
: Cit ASPUGAOD 
re = 
de ‘ith '- of FY 
’ JL : - >“ 


* oa 
FERIBG: Siscupe @0f67erb pabbeoaee 


i ‘ 7 


4 oe Se 7 che £63 22 f evip ot 


_ 


+ Ad j me ia SF. aye a 7 y 
- 4 
-? et 5 ae a 
vs rits < De L1O0 fIG ai 7 oO SG oh ID Ao Lae 1 “ae 
’ 7 <: 7 : a 


ASV od inew eb 44 


30 


multiple coverage of the same subsurface points with 
different shot and detector positions. The effectiveness of 
this method is dependent on applying the proper time 
corrections to each trace so that the primary reflections 
wili be moved in phase and properly stacked. One reguires, 
as a result, a method to calculate corrections for all the 
traces with their different shot point to geophone 


distances. 


For a Simple horizontally layered model these time 


corrections could be computed from the equation 


T(x,n)2 = T(o,n)2 + X2/V(n)2 (355) 


T(x,n) is the two-way travel time at offset X for 


layer n 
T(o,n) ig the two-way normal incidence time for 
layer n 
V(n) is the RMS velocity as given by Dix (1955) 
V(n)2 = SUM{v (i) 2t (i) } /SUM{t (i) } (3.2) 
SUM indicates summation over index i from 1 ton 
v(i) is interval velocity of the i'th layer 
t(i) is two-way normal incidence time in i'th 
layer 


Now, the two-way normal incidence time for layer n is 


T(O,0) = SUM {t (2)} (3.3) 
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and the interval velocity for layer n is 


WO eee (ty) St Oe tf) en (oO; n= 1) | Aceh (oO, 0) =f (0,n- 1) ) 


(3.4) 


Rederivationsof cequationss 3.41° and)i3.2 is ‘given -in 
Schneider and Backus (1968) and Taner and Koehler (1969). 
in the presence of dipping beds, model experiments have 
Shown that the time-distance relations still remain nearly 
hyperbolic. However, the RMS velocity now obtained is only 
an apparent stacking velocity and should not be used to 
calculate interval velocities unless corrections are made 


for dip (Taner et al, 1970). 


In general, the method of velocity determination from 
refiection data is done on a grid system. A value of T(0o,n) 
is chosen and V(n) iS varied between some minimum and 
Maximum. Each value of T(o,n) and V({n) defines a particular 
set of time corrections. These corrections are applied to 
the data and a coherency measure is computed for each time 
correction. The final display shows this coherency measure 
plotted as a function of incidence time and RMS velocity. 
Peaks in the display indicate arrivals cf coherent energy at 
particular values of incidence time and RMS velocity. Such 


a display is referred to as a velocity spectrum. 


The coherency measure is computed within some time gate 


centred about the time corrections for each trace as 
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Calculated by equation 3.1. It may take many forms from 
Simple summing or normalized summing of traces to more 
complicated cross correlation techniques. Garotta and 
Michon (1967) use summing of traces as a coherency measure 
while Schneider and Backus (1968) describe a technique based 
on cross correlations. Taner and Koehler (1969) give two 
coherency measures. One is a sum of cross correlations 
between pairs of traces while the other is a function which 
they cali the semblance function. A gocd review of various 


coherency measures is given by Montalbetti (1971). 


Robinson (1969) describes a high resolution velocity 
technique which performs the calculation of a velocity 
spectrum in the frequency-wavenumber domain. in ccmparison 
to a cross correlation method in the space-time domain this 


method is faster On a computer. 


3.2 The Velocity Pregram Used in this Thesis 


If the common depth point method is not used for 
recording the data, the velocity spectral methods discussed 
in the preceding section are still applicable for uniforn, 
flat layering. However, if significant dips are present on 
the interfaces, the velocity spectra can not be used for 
determination of interval velocities. In such a case it is 
necessary to find a different approach to velocity 
determination using the computer. Due to the severe 


restrictions imposed on shotpoint and detecter locations by 
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the extensive farming activity in the area it was not 
possible to arrange the profile to conform to the usual 
requirements of the common depth method. Further, a 
preliminary analysis of the data indicated that large dips 
were present on the reflecting horizons. Therefore a 
velocity program was written specifically for application to 
the particular set of data that was used in this thesis. 
This program was designed for determining velocities, depths 
and dips of layers within the basement and requires a_ known 
velocity model for the region above the basement. Models 
consisting of one or two layers within the basement may be 
treated although the program could easily be generalized to 
handle more layers. Straight-line ray paths within the 


layers have been assumed. 


The program can be briefly summarized as follows. A 
particular velocity, normal incidence time and dip wiil 
specify one model of the crust. Knowing the shot to 
detector distance one can solve iteratively for the ray path 
corresponding to this distance for this particular model. 
This allows calculation of the travel time for each shot to 
detector distance. A time correction is applied to each 
trace on one record and a measurement made of the coherency 
of the energy arriving along this time correction. By 
systematically varying the velocity, incidence time and dip 
one calculates the coherency measure as a function of these 
three parameters and peaks in this display will indicate in- 


phase arrivals of energy. 
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The detailed operation of this program can be 
understood from examination of figure 3.1 which illustrates 
the case where there is only one layer below the basement. 
The ray path is ABCDE. The sections AB and DE represent the 
path through the layers above the basement. The parameters 
which are specified are the velocity in the first layer, Vi, 
the dip of the boundary at the bottom of the layer, c1, and 
the normal incidence time for this layer, 11. The normal 
incidence time is the time it would take a P wave to travel 
from a point directly below the shot to the bottom of the 
layer and back along a path perpendicular to the bottom of 


the layer. 


The program assumes that the top of the basement is 
flat which is consistent with the velocity model and 
formation depths frcem wells in the area as described in 
Chapter 2. Using the velocity modei for the above-basement 
region and assuming straight line paths in the various 
layers of that model one can, with the aid of Snell's law, 
calculate travel time and x-distance travelled by a ray 
making an angle a1 with the normal just above the basement. 
Table 3.1 gives travel times and x-distances for various 
values of a1. For the purpose of the program, values for x- 
distance are stored in a table for integral values of al and 
interpolation is used to obtain an exact value when ail is 
not an integer. Sufficient time accuracy for any angle 


between two integral values of a1 may be obtained by using 
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Figure 3.1 


One layer model diagram. 
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TABLE 3.1 


Travel Time (sec) 


0.724 
0.724 
0.724 
0.725 
Osa 
0.726 
0.727 
0.728 
0.729 
0.730 
0.731 
0.733 
0.734 
0.736 
0.738 
0.740 
0.742 
0.744 
0.747 
0.749 
aloe 
0.755 
0.758 
0.762 
0.765 
0.769 
0.773 
0.777 
0.782 
0.786 
054% 
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X-distance (meters) 


628 
667 


a2 
954 
S93 
1042 
1088 
1134 
1181 
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the value of time for an angle half way between the integral 
values. The tables of x-distance and time as a function of 
al are referred to as UPX(a1) and UPT(a1l). 


In reference to figure 3.1 


BY = PQ + BPeSIN(c1) (3.5) 
BC = BY/COS (a2+c1) (3.6) 
CM = BC*COS (a2) (3.7) 
BM = BCeSIN (a2) (3.8) 
MD = CMeTAN(a2+2c1) (3.9) 
CD = CM/COS (a2+2c1) (3.10) 


The x-distance travelled in the below-basement layer is 


X' = BD = BM + MD (3.11) 


PLOMmeS.0 sol, 2.0, and 3.9 


BD = BY[ {SIN(a2) +COS (a2) #TAN (a2+2c1)}/CCS (a2+c1) J] 
(3.12) 


The travel time below the basement is 


g' = BC/V1 + CD/V1 (3.13) 
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Pron 86<.6,0e% 7Vand?sS.10 


T' = BY/V1[ {1+COS (a2) /COS (a2+2c1) } /COS (a2+c1) } 
(3.14) 


if we definesythe fLunctaons*Fxri(ci;a2) and FT(cl,a2)+ to 
be the parts of equations 3.12 and 3.14 respectively within 


the curled brackets, { }, then 


X' = BYeFX(c1,a2) (3.15) 


T!' = (BY/V) eFT(c1,a2) (3.16) 


Values of the functions FX and FI can be tabulated for 
integral values of cl and a2. For values of these angles in 
the range of interest, linear interpolation gives sufficient 
accuracy (1 msec) for non-integral values of the angles. 
Appendix A gives some sample values cf these functions. 
Now, if we specify the velocity, V1, and incidence time, T1, 


then 


PQ = VieT1/2 (3217) 
in addition, 

BP = UPX (a1) (3.18) 

DT =sUPX (b1) (3.19) 
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But we know the velocity in the lowest layer above the 
basement, VL, and from Snell's law we can calculate a1 and 


bi. 


a1 = ARCSIN[ VLeSIN(a2)/V1] (3.20) 


b1 = ARCSIN[ VLeSIN(a2+2c1)/V1] (3.21) 


The total x-distance and travel time are 


to Li Bp DT (3222) 


T = T' + UPT(a1) + UPT(b1) (3.23) 


It is apparent from equations 3.5 and 3.15 to 3.23 that 
if Vi, 11 and cl are specified then X and T depend only on 
a2. The progragr starts by calculating X for an estimated 
a2. If this X doesn't match the true shot-recorder 
distance, within a suitable distance (25 meters), then a new 
estimate of a2 is calculated based on the previous estimate 
andes thes error in. -X. After this iterative procedure has 
yielded the correct a2, the travel time for this path is 
calculated. This is done for all shot-recorder distances on 
a particular record. These time corrections are then used 


in calculating a coherency measure. 


The coherency measure is calculated for a time gate 


centred about the time correction. Time gate width is twice 
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the value of the increment on incidence time. The coherency 


measure used here is 


C = SUMt [SUMi ff (i,t)} J2 
SUMt [SUMi{{f (i,t) {} ]2 (3.24) 
where SUMt is summation over t 
SUMi is summation over i 
i is the index on the number of traces 
t is the time index within the time gate 


£(i,t) represents the seismic trace 


The numerator of equation 3.24 is a measure of the 
stacked energy within the time gate and the denominator is a 
measure of the vaiue that the numerator would have if all 
traces being stacked were in-phase. C approaches zero for 


out of phase arrivals and one for in-phase arrivals. 


By systematically varying V1, T1 and cl one caiculates 
the coherency measure as a function of the three variables 
and maxima indicate arrivals at particular values of V1, fT1 


and) cit: 


For the two layer case the argument is the same as for 
the one layer case except that now V1, 11 and ci are known 
quantities from the previous calculation and the 
corresponding parameters V2, T2 and c2 in the second layer 
are varied. In a similar manner the argument could be 
extended to three layers. Appendix A ccntains the equations 


for the two-layer case and Appendix B has a listing of the 
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CHAPTER 4 


INTERPRETATION 


mt tntroductidn 


Ali of the data was plotted and examined for good 
reflections. The plotted records were usually filtered with 
a5 to 30 Hz four pole Butterworth band-pass recursion 
filter (Shanks, 1967). As an example of a filtered record, 
figure 4.1 illustrates shot 1 from September 2. Each trace 
on this plot is normalized to the same maximum for plotting 
Simplicity and the two segments have been normalized 


separately. 


The velocity program used to analyze the data has been 
described in chapter 3, and the velocity modei for the 
sedimentary section has been described in chapter 2. From 
figure 4.1 we can see that there are two prominent 
reflections at times of 7.5 and 10.7 seconds and a third 
weaker reflection at 11.7 seconds. These prominent 
reflections are common to most records and thus they have 
been chosen as the two layers to be used in the velocity 
model. An attempt was made to include a layer corresponding 
to the reflection at 11.7 seconds, but because of the poor 
resolution due to the iarge depth and the small horizontal 
displacements, the velocity analysis was inconclusive for 


this reflector. Although a velocity could not be obtained 
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for the 11.7 second reflection, one can conclude that the 
deepest layer must range in thickness from 3 to 4 km for any 


reasonable velocity. 


The records were checked for good reflections in the 
first 1.5 seconds, corresponding to the sedimentary section, 
so that a possible comparison could be made with a synthetic 
seismogran. Unfortunately, the gain settings used in 
recording the data were such that this section of the record 


was unusable. 


Because of the suggestion of a low velocity layer in 
the upper basement by Sprenke (1972), it was felt that a 
careful analysis of the data in the appropriate time range 
should be made. A velocity spectra program as described by 
Taner and Koehler (1969) was applied to the data from the 
shots with small horizontal displacements. This program 
assumed horizontal layering. The results of this study 
indicated that prominent reflections to 3.5 seconds have RMS 
velocities in the range of 2.8 to 3.5 km/sec which are 
consistent with multiples from the sedimentary section. 
Using the velocity model for the sediments, a synthetic 
seismogram (Clowes, 1969) was produced. This synthetic 
seismogram showed multiples to 3.5 seconds and matched the 
field records reasonably well. For these two reasons it was 
felt that there is no justification in the present data for 


any suggested interfaces in this region of the crust. 


For the above reasons, it was decided to construct a 


; ; Bs 
iz Jn244 HU Lona 3 i s. m0 ,7oLtoe lie moll 


+ ee o a 
Jo GA fF - & 3 2090290587 BRE SVUBI ava daeuk Faeq4s ou “a 
t 7 
a 
ptiooisv, eiltsacaset : 
; 108 hedaw 19 2brooe2 dy == a 
; } rl OD ,80gooe@e €.f Faxkt 
: » il id 3 f }flaseou ¢ saat oe 
r . 49 os 4 «4 4 Oe ropa sili a bers 
f . 7 te Hh of ‘rit = el #3 i i> limes p1008s | 
n lds 
“ F 
6 i ae t t 13 20 
ban 
2 ae Pr) } d Yi 38984 z 
3.2 : ras 3 €zavipas 
~— 
i I r Fi0GCLSY A +855 4 


~ ‘ > { * - 
4 a - & £ eee | 1 be he OS ‘sh 
- r 7 
z 3 ovn.sdeysi Lstaqosiuod 
4 - —_ al 
. > » abe 4 ¢ a q tern ne s &» mn > wg | 
7° sav lyos 1 2890 1GGIG 7447 Sar eQEoRe 
: its 


+r. i of Bes de spit eit ak sektioe 
/ a ab 


* 
- 


oli 43 ~ . : ; : 14 7 < atch Fs —p | ;isiea- n+. $2 te Lee 


siaomihas® any Goi Lebos yehoolgy 2a) pt 20 
aa2e8 G ew 


a 
o 


; ; 
voge 7 nstz0 Ss98 


45 


two layer model of the crust. As an example of the output 
from the velocity program, table 4.1 and figure 4.2 are 
included. The table contains results obtained by running 
the velocity program on the 7.5 second reflection for shot 2 
on September 2. These same results are Shown as contours in 
the diagram in figure 4.2. The velocities are in km/sec and 
the times are two-way normal incidence times in seconds. 
The numbers in the table range from zero to one thousand 
depending on the coherency of the arrivals along various 
time corrections. From these results it can be seen that 
dip and velocity can not be completely separated, since 
increased dip results in a shift of the peak to higher 
velocity. To help improve on the interpretation, calculated 
travel times were printed by the velocity program for the 
various models. These are compared to times read from the 
records and this aids in the choice of a final model. For 
the example shown, the velocity and dip picks would be at 


6.6 and 15 respectively. 


4.2 The Iwo Layer Velocity Model 


The September 2 data consisted of two groups of three 
shots at opposite ends of the spread. These shots were 
independently interpreted and the final model was based on 
the average of velocity and dip for each group. The three 
shots to the north gave a layer which dipped at 15 degrees 


to the south with velocity of 6.4 km/sec and two-way normal 
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incidence time in this layer below shotpoint 2 of 5.8 
seconds. The shots to the south led to a layer which dipped 
to the south at 18 degrees with a veliocity of 6.3 km/sec and 
incidence time below shotpoint 5 of 5.69 seconds. No 
attempt has been made to phase correct the times to the 
beginning of the reflection. Estimated uncertainties are 
less than 2 degrees in dip, .2 kmysec in velocity, and .05 
seconds in time. These uncertainties are consistent with 
scatter in individual determinations of dip and velocity. 
These two reflecting layers are shown in figure 4.3. The 
normal incidence times have been converted to depths and the 
horizons are plotted at the appropriate subsurface 
positions. It is obvious that the reflections are either 
from different discontinuities or that a fault exists in the 


reflecting interface. 


Although velocity resolution was poor on the east-west 
profiie, it was possible to choose a reasonable velocity 
from the above results. Dip must then be between 10 and 15 
degrees to the east. This dip would introduce some error in 
interpreting the structure along the north-south profile 
because the velocity program makes the assumption that no 
dip is present in the east-west direction. It was felt that 
the data was not adeguate enough to justify a correction for 


this effect. 


Based on the resuits discussed in the previous 


paragraph, the second layer was added. The shots to the 


ry - 
sabao | y= z 
\ t ® ¢fioolsy 4 d4#iw @esapeh Sf ga datce odt oF 
, ws 1 ee 

: . : ; at 4 | 
> ; © tihogtede woled esis ssasbingt 
aan” ° 
3 ig of shea feed eed teaes98 


peer fF Lal Perk otras Lie2z ef? > paigatosd 
j sent seal 


T .sdRit AE BORODaS 
7 _. 


* ‘ ; 5 
f oO 2 GiLle Feo isnebivibas lt : 
f i oa I 3S 4 BHlsoosgise ws = 
+ + 
~~ ' ‘ 
a ‘ r ¥v aye r enn Wolds 
~~ 


- ; ' e - * po ~ P a . 4 J - 

ae a ; bayoaggs oid erztole: xs 21OS L108 
a”. 
2 . a a a ina. ~4 . anos: +?teee 
= £FDS4 i. 25) x VOLVEO B22 Jz »-2n0z£tLaC 7d 
aan 


i 


- as ; x ra 

ted3 10 eeitiupsiaopebh. Ynaoasitlb soa 
2 ha 

,sostistut palsoslieg 


~ Siz it ao rocg 264 S0fsuloaes yetogiew desalsia 65 


i) 


. J eidiarecy sev. 2s! 
i ets 47 Jens gid er lees evods edt BORE 


2 


101% Du iyow Gin eBLaT s2259 S23 OF asgzpob 


= 
4 


& 
Semg dfpoe-d320n ei 3 Wess S7UIOUBZIE Sis wah 
1 ae 7 


ISAT Gosigaveas ods Sorts MSIp poz, el ioete oat ast 
ae ol r S90 
anv -- +M9is08zib Jeev-se 282 ad ‘at sa el i 
oe ey “9 
= ~ al if Ja : ‘ ; 
. 7 are Ss «&  V¥, bts uC a3 ‘¢ one — ‘ie wy ; “2 ase ; 


pi, _ _ 


e 
{> 4-3 


hi 


south of the spread showed no clear reflection at the 
appropriate times, and hence the second layer is based only 
on the three shots to the north. The model obtained for 
this second layer has a velocity of 6.5 km/sec, a dip of 2 
degrees to the south and an incidence time of 4.15 seconds 
in this layer under shotpoint 2. The uncertainties in these 
values would be slightly greater than in the values for the 
first layer because this layer depends on the previous 


results. This second layer is also shown in figure 4.3. 


The above reflection was absent on the 7 traces at the 
north end of the recording spread for shot 3 and on the 3 
traces at the north end for shot 2. This could be explained 
by the presence of a fault in the upper layer. fhe ray 
paths for the reflections which are absent would pass 
through the region of the fault and be obscured. This would 
substantiate the previous suggestion of a fault in this 


region. 


4.3 Conclusion 


The crustal model found for this area has shown the 
presence of steeply dipping interfaces within the deep crust 
and has suggested a possible fault. It would be desirable 
to have a continuous profile across the location of the 
suggested fault and greater offset distances would be 
advantageous for the deeper sections. Unfortunately, in 


this’ “region ‘the’ presence *‘ofs pipe lines and power lines 
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Figure 4.3 Crustal Model Diagran. 
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severely limits the areas of study. One solution to this 
problem might be the use of vibroseis methods; however, this 
might not have sufficient depth penetration to be useful. 
The presence of dipping layers within the basement suggests 
that P-coda results which were based on models consisting of 


flat crustal layers should be re-examined. 


It is interesting to compare these results with those 
obtained in southern Alberta by Kanasewich and Cumming 
(1965), Clowes et ai (1968) and Clowes and Kanasewich 
(1970). Kanasewich and Cumming (1965) mention a weak 
reflection at 8.4 seconds which is certainly less prominent 
than the 7.5 second reflection here. The 10.7 second 
reflection could be identified as the R discontinuity on the 
basis of its general appearance. This reflection is .7 
seconds eariier than in the south which is consistent with 
the fact that the 7.5 second reflection is also earlier than 
the corresponding reflection in the south. If) sthei etis7 
second reflection is interpreted as the Moho then the 
interval thickness from the R to M discontinuities is much 
less than in the southern part of the province. The total 
crustal thickness wouid then be 35.5 km, assuming a velocity 
in the lowest layer of 7.2 km/sec as found in the southern 
part of the province. This crustal thickness is less than 
in southern Alberta and the main reason would be the 
apparent thinning of the deepest layer between the R and M&M 


discontinuities. 
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APPENDIX A 


TWO LAYER VELOCITY PROGRAM 


The derivation of the equations for the two layer 
velocity program is very similar to that for the one layer 
velocity program described in Chapter 3. In figure A.1 the 
ray path is ABCDEFG. PQ is normal to the first reflecting 
surface and is specified by V1 and 11. RS is normal to the 
second reflecting surface and R is directly below the shot 
point. V2, and T2 specify RS. The angle a3 is the variable 
used in iteratively solving for the ray path and ali other 
angles can be found from it. Referring to figure A.1 and 
letting VL be the velocity in the lowest layer above the 


basement we have: 


a2 = ARCSIN{V1¢eSIN (a3) /V¥2} - c1 (A. 1) 

al = ARCSIN{VL®SIN (a2) /V1} (A. 2) 

b2 = ARCSIN{V1eSIN(a3+2(c2-c1))/V2} + cl (has) 

b1 = ARCSIN{VLeSIN(b2) /V1} (A.4) 

Now, BP = UPX (a1) (A. 5) 
FK = UPX(b1) (A.6) 
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Figure A.1 


{Two layer model diagraa. 
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CZ 


CD 


EM 


DE 


PrOMW Aw IS s. be lo, 


CE = CZ[SIN(a3)_+_COS (a3) eTAN (a3+2(c2-c1)) J 


PQ = V1eT1/2 


= PQ + BPeSIN(c1) 


= BY/COS (a2+c1) 


BV = BCeSIN(a2) 


= {BP+BV}/COS(c1) 


RS = V2eT2/2 


RS + CReSIN(c2-c1) 


=) CL/COS(a3tCc2—-c1) 


DM CDeCOS (a3) 


CM 


CDeSIN (a3) 


DMe TAN (a3+2 (c2-c1)) 


DM/COS (a2+2 (c2-c1) ) 


A.16 and A.1/ 


COS fastc2-c1) 
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(A.7} 


(A. 8) 


(A.9) 


(A. 10) 


(A. 11) 


(A.12) 


(A. 13) 


(A. 14) 
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Also, WV = CEeCOS (c1) (A.20) 


The travel time in the bottom layer is 


T" = CD/V2 + DE/V2 (A. 21) 


PEOM A StG, €f..1 5hand $A ats 


m= C€2(1_+ COS (a3) /COS (a3+2 (c2-c1)) J 
V2eCOS (a3+c2-c1) (A.22) 


Recalling the definitions of FX and FT in 3.12 and 3.14 


CE = CZeFX(c2-cl,a3) (A.23) 

Tt = | -CZ/V 2 jeFT (c2-¢1,a3) (A. 24) 

Now, EW = PQ/COS(c1) + (BP+BV+WV) eTAN (c1) (Ae) 
EF = EW/COS(b2) (A. 26) 

FW = EWeTAN(b2) (Av27)} 


But, the x-distance is 


xX =. BP + BY 4°00 + fe + FR (A. 28) 
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and the travel time is 


T = UPT(a1) + UPT(b1) + (BC+EF)/V1 + T" (A. 29) 


Prom equarions A. hitOen slay Awe, (Aveo, Acoo, As2i ‘and 
A.28 it is clear that when V2, T2 and c2 are fixed along 
wicnroavi, 01 angecl;ethen * 15a function of only a3. For any 
value of X it is possible to calculate the travel time for 
the corresponding path according to eguations A.1 to A.13, 
A.20, A.23 to A.Z6 and A.29. Thus, in the two layer case one 
can solve iteratively for the value of a3 which gives the 
proper x-distance. This can be done for all shot-recorder 
distances on a record and time corrections can be calculated 


for each trace in the same way as for the one layer case. 


Table A.1 which follows gives values of the functions 


FX and FI for values of c2-c1 within the range of interest. 
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sae 
26.0 
27.0 
28.0 
29.0 
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34.9 
32.0 
33.0 
34.0 
35.0 
36.0 
S720 
38.0 
39.0 
40.0 
471.0 
42.0 
43.0 
44.0 
45.0 
46.0 
47.0 
48.0 
49.0 
50.0 


FX 


0.90 

0.0369 
0.0734 
0.1095 
0.1451 
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0.2155 
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0.2846 
e318? 
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FX 


0.0 

0.0349 
0.0698 
0.1048 
0.1399 
0.1750 
0.2102 
0.2456 
0.2811 
0.3168 
0.5527 
0.3888 
0.4251 
0.4617 
0.4987 
0.9359 
0.5735 
0.6115 
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0.7279 
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2.0000 
2.0003 
2-0012 
2.0027 
2.0049 
2.0076 
2.0110 
20050 
2.0197 
2.0249 
2.0309 
2.0374 
2.0447 
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0.6563 
0.7000 
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2.0065 
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2.0214 
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2.0490 
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2.1560 
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2.3485 
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2.4045 
224346 
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2.4995 
2.5344 
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1.1309 
Tego 
1.2621 
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APPENDIX B 


LISTING OF VELOCITY PROGRAM 


This appendix lists the actual velocity program. Table 
B.1 lists the correspondence between variables used in this 


thesis and the variabie name in the progran. 


The program is in Fortran and for maximum efficiency 
should be compiled using the IBM Fortran H compiler with the 
optimization level set at 2. A typical run for 12 seismic 
traces with 20 different velocities, 10 different incidence 
times, 5 different dip angles and a time gate width of .1 
seconds required approximately 15 seconds on an IBM 360/67 


computer. 
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THESIS VARIABLE 


UP® 


UPX 


FX 


Fz 


at 


b1 


a2 


a3 


b2 


pe 


cz 


v1 


v2 


T1 


T2 


TABLE B.1 
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PROGRAM VARIABLE NAME 


UPFUNT 
UPFUNX 
FUNAD 

FUNAED 


ALPHA1 


ALPHA (1,1) 
ALPHA (2,1) 
BETA (1,1) 


THETA1 or 
THETA for 1 layer 


THETA (2 layer case) 


V1 or 
V for 1 layer case 


V for 2 layer case 


ToL 
TO for 1 layer case 


TO for 2 layer case 


CARR CI CK CR CII OR AC GOR OI I I i a ak ak a ak ak ak ak ak ak ak ak ak 


CER EKER EE TWO LAYER VELOCITY PROGRAM FEET E ERK EE 
CR RR aI ICCA ACI II I a ROK a ak ak 2k 3k 2k ak ak ak 


THEY PURPOSE? OFS THIS) PROGRAM® ISVTOUTESTI ASSET OF 12 
SEISMIC TRACES FOR AN ARRIVAL OF A SIGNAL WITH A 
PARTICULAR STEPOUT IN TIME FROM TRACE TO TRACE. THE 
FACTORS AFFECTING STEPOUT ARE DIP VELOCITY, AND 
DEPTH OF REFLECTOR. 


THE PROGRAM IS DESIGNED UPON THE ASSUMPTION THAT A 
KNOWN FLAT LAYERED MODEL EXISTS ABOVE THE BASEMENT AND 
THAT WE WILL LOCK FOR NO MORE THAN 2 LAYERS BELOW THE 
BASEMENT. EACH OF THESE 2 LAYERS IS FOUND BY A SEPERATE 
RUN OF THE PROGRAM STARTING WITH THE LAYER JUST BELOW 
THE BASEMENT AND GOING DEEPER ON EACH SUCCESSIVE RUN. 
WHEN WORKING ON THE SECOND LAYER BELOW THE BASEMENT THE 
RESULTS OF THE FIRST LAYER ARE NEEDED. THE QUANTITIES 
V, THETA, AND TO IN THE UNKNOWN LAYER ARE VARIED AND 
TIME CORRECTIONS ARE CALCULATED FOR EACH TRACE FOR EACH 
DIFFERENT SET OF V, THETA, AND TO. 


V = INTERVAL VELOCITY IN LAST LAYER 

THETA = DIP OF LAST LAYER 

TO = TWO WAY TRAVEL TIME ALONG A PATH IN THE 
LAST LAYER PERPENDICULAR TO THE INTER- 
FACE AT THE BOTTOM OF THIS LAYER AND 
FROM A POINT DIRECTLY BELOW THE SHOT. 


THE ARRAYS UPFUN, UPFUNT, FUNAD, AND FUNAED ARE 
REQUIRED. THESE ARE USUALLY STORED ON A DISK (IN THIS 
CASE LOGICAL UNIf£ 1). 


UPFUN (61) CONTAINS THE X-DISTANCE (DISTANCE ALONG 
GROUND FROM SHOT POINT) WHICH A SEISMIC RAY TRAVELS FOR 
A GIVEN ANGLE BETWEEN THE RAY AND THE PERPENDICULAR TO 
THE INTERFACE IN THE LAYER ABOVE BASEMENT. THE ANGLE 
CAN RANGE FROM -30 TO 30 DEGREES AND X VALUES ARE GIVEN 
IN UPFUN FOR ALL INTEGER VALUES IN THIS RANGE. EXACT 
VALUES ARE FOUND BY INTERPOLATICN. 


UPFUNT (61) IS SIMILAR TO UPFUN ONLY IT CONTAINS 
THE TRAVEL TIMES TAKEN BY THE RAY ABOVE THE BASEMENT. 
IN THIS CASE INTERPOLATION IS NOT USED BECAUSE 
TRUNCATION IS ACCURATE ENOUGH. UPFUN(1) IS THE TIME FOR 
AN ANGLE OF -— 30.5, UPFUN(2) FOR. -29.5 ... UPFUN (31) 
FOR -.5.0R .5 AS BOTH ARE EQUAL, UPFUN(32) FOR 1.5 ETC. 


Seer or@e One GOO Oi AAA AA AAA AAR ARO AMADA AAAAAMA AQAA AaAaAan aan @Qana 
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FUNAD (I,J) CONTAINS VALUES SUCH THAT FUNAD (I,J) 
MULTIPLIED BY THE THICKNESS (PERPENDICULAR TO BOTTOM OF 
LAST LAYER) OF LAST LAYER AT THE POINT THE RAY ENTERS 
THE LAST LAYER EQUALS THE X-DISTANCE TRAVELLED IN THE 
BOTTOM LAYER. THE I INDEX IS THE INDEX ON THE DIP 
ANGLE. IF THETA IS THE DIP OF THE BOTTOM LAYER AND 
THETA] IS THE DIP OF THE LAYER ABOVE IT THEN I=1 
REPRESENTS (THETA-THETA1) EQUAL TO -20 AND I=2 IS 
(THETA-THETA1) EQUAL -19 ETC. THETA AND THETA1 MUST 
TAKE ON INTEGRAL VALUES. THE J INDEX IS THE INDEX ON 
THE ANGLE THAT THE RAY MAKES WITH THE NORMAL WHEN IT 
ENTERS THE BOTTOM LAYER. J=1 IS A VALUE OF THIS ANGLE 
EQUAL TO - (THETA-THETA1) AND J=2 IS 1- (THETA-THETA1) 
ETC. INTERPOLATION CAN BE USED. NOTE THAT IF SHOTS ARE 
SHOT IN AN UPDIP DIRECTION THEN THETA IS NEGATIVE. 


FUNAED (I,J) IS SIMILAR TO FUNAD ONLY FUNAED (I,J) 
MULTIPLIED BY 1/V TIMES THE MULTIFLYING FACTOR FOR 
FUNAD(I,J) GIVES THE TRAVEL TIME IN THE BOTTOM LAYER. 


USING THESE ARRAYS AND A GIVEN V, THETA, AND TO IT 
IS POSSIBLE TO CALCULATE A TIME CORRECTION FOR EACH 
TRACE WITH ITS PARTICULAR X-DISTANCE. 


V, THETA, AND TO ARE ALL VARIED IN AN ATTEMPT TO 
FIND THE VALUES WHICH BEST FIT THE PARTICULAR 
REFLECTION ON THE RECORD. AS A MEASURE OF GOODNESS OF 
THE FIT THE FOLLOWING IS CALCULATED: 


SUM(T) = SUM OF F(T) OVER ALL TRACES INSIDE 
A TIME WINDOW CENTRED ABOUT THE TIME 
CORRECTION. THE WINDOW WIDTH IS 
TWICE THE INCREMENT ON TO. 


SUMNOR(T) = SAME AS ABOVE ONLY ONE USES 
ABSOLUTE VALUE OF F(T). F(T) IS 
THE SEISMIC SIGNAL AS A FUNCTION 
OF TIME. 


THEN SUMX = SUM OVER T OF SUM (T) **2 
SUMNX = SUM OVER T OF SUMNOR (T) **2 


SUMX IS A MEASURE OF THE TOTAL STACKED ENERGY 
COMING THROUGH THE WINDOW. SUMNX IS A MEASURE OF WHAT 
THE TOTAL STACKED ENERGY WOULD BE IF ALL TRACES BEING 
STACKED WERE PERFECTLY IN-PHASE. THE QUANTITY 
1000*SUMX/SUMNX WILL APPROACH ZERO FOR RANDOM ARRIVALS 
AND WILL APPROACH 1000 IF ALL ARRIVALS ALONG THE TIME 
CORRECTION ARE IN PHASE. 
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Cc 
DIMENSION UPFUN (61) ,UPFUNT (61) ,FUNAD (45,31) , FUNAED (45, 31) 
DIMENSION A(500,12),X(12) ,TCOR(12) , ALPHA (2,12) ,BETA (2,12) 
DIMENSION ISPEC (20) ,SUM(40) ,SUMNOR(40) , TITLE (20) , VEL (20) 
CK KK 
fe TITLE = TITLE OF THIS RUN 
Cc NUMTR = NUMBER OF TRACES FOR THIS RUN 
G IP = POINT IN ARRAY ON TAPE THAT WE WANT AS FIRST POINT 
€ OF 500 FOR SEISMIC TRACES STORED IN A 
Cc LAYNUM = NUMBER OF LAYER (BELOW BASEMENT) BEING SOLVED 
Cc V1 = VELOCITY IN FIRST LAYER BELOW BASEMENT 
c T1 = TO VALUE FOUND FOR FIRST LAYER BELOW BASEMENT 
Cc ITHET1 = DIP OF FIRST LAYER BELOW BASEMENT 
Cc X = ARRAY OF X-DISTANCES OF GEOPHONES 
c TINIT = INITIAL TO VALUE 
Cc TIMINC = INCREMENT ON TO 
G NTIME = NUMBER OF TO VALUES TO BE USED 
fe VINIT = INITIAL VELOCITY VALUE 
(S VELINC = VELOCITY INCREMENT 
Cc NVEL = NUMBER OF VELOCITY VALUES TO USE 
Cc ITHETA = INITIAL VALUE OF DIP ANGLE 
eC ITHINC = INCREMENT ON THETA 
Cc NTHETA = NUMBER OF THETA VALUES TO USE 
C** KX 
READ (5,1) TITLE 
READ (5,2) NUMTR,IP,LAYNUM,V1,11,1THET1 
READ (5,3) X 
READ (5,4) TINIT, TIMINC,NTIME, VINIT, VELINC,NVEL,ITHETA 
e, 1THINC, NTEETA 
NGATE=357. 1428*TIMINC+.5 
IF (NGATE.GT.40) STOP 1 
IF (NVEL.GI.20) STOP 2 
IF (LAYNUM.GT.2) STOP 3 
WRITE (6,8) TITLE, (X(I) ,1=1,NUMTR) 
WRITE (6,9) LAYNUM 
GC TO (44,43) ,LAYNUM 
43 WRITE (6,10) V1,T1,ITHET1 
THETA1=ITHET 1/57. 29578 
STH1=SIN (THETA1) 
CTH1=SQRT (1. -STH1**2) 
VT21=V1*T1/2.0 
GO TO 45 
44 LTTHET1=0 
45 VEL (1)=VINIT 
DO 50 I=2,NVEL 
50 VEL (I) =VEL (I-1) +VELINC 
READ (1) UPFUN,UPFUNT,FUNAD,FUNAED 
CK KK 
Cc RGET READS IN SEISMIC TRACES OFF TAPE INTO A(500,12) 
Cc WITH A(1,I) EQUAL TO THE IP'TH POINT OF THE I'TH TRACE 
Cc WHERE THE TRACE ON TAPE IS ASSUMED TO HAVE POINT 1 AT 
Cc TIME ZERO 
C** KX 


CALL RGET(2,11208,A,500, NUMTR,IP,0,&999) 
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CK KK 
THREE LOOPS FOLLOW: 
FCT IIIS EOOPSIN PWHIGH THETA IS VARIED 
LET2Z 01S LOOP IN WHICH TO IS VARIED 
terse 0S LOOP EN WHICH “¥ -IS* VARLED 


Cc 
Cc 
Cc 
Cc 
¢ 
Cc ALPHA REFERS TO ANGLE BETWEEN RAY AND NORMAL FOR 
¢ DOWNGCING RAYS. 
Cc ALPHA(I,J) = ANGLE AT TOP OF I'TH LAYER BELOW BASEMENT 
Cc FOR J'TH TRACE 
Cc BETA IS SAME AS ALPHA ONLY FOR UPGOING RAYS 
Cc ALPHA1 AND EBETA1 ARE ANGLES IN LAYER ABOVE BASEMENT 
Cc XP = X DISTANCE TRAVELLED BY RAY BEING CONSIDERED 
C**K KX 
DO 500 IcT1=1,NTHETA 
THETA=ITHETA/57.29578 
STH=SIN (THETA) 
CTH=SQRT (1.-STH**2) 
STHM1=SIN (THETA-THETA 1) 
I1=ITHETA-ITHET 1421 
TE (2176T24520R. 11.1821) STOP 4 
TO=TINIT 
WRITE (6,5) TITLE, ITHETA 
WRITE (6,6) (VEL(ZI) ,I1=1,NVEL) 
DO 400 ICT2=1,NTIME 
V=VINIT 
DO 300 ICT3=1,NVEL 
VI2=V*T0/2.0 
CK KK 


Lé: COMPUTED GO TO BRANCHES TO 2 PLACES DEPENDING CN NUMBER 
€ OF LAYERS INVOLVED. THESE 2 CASES WERE CODED SEPARATELY 
Le AND THEREFORE ARE NOT COMBINED AT ALL. THEY PROBABLY 
C COULD BE CCMBINED TO SOME EXTENT TO REDUCE CODING 
Lé AND STORAGE. TO SAVE STORAGE ONE CAN JUST OMIT THE 
(a: TWC UNUSED SECTIONS WHEN COMPILING. 
Cc 
Cc DUE TO A BUG IN THE FORTRAN H COMPILER THESE 2 SECTIONS 
Cc CAN NOT BE COMPILED TOGETHER. TO COMPILE ONE LAYER AT 
Cc A TIME REPLACE STATEMENT 110 OK 140 BY STOP 110 OR STOP 
G 140 STATEMENTS AND LEAVE OUT THE ENTIRE SECTION. 
C#K KK 
GO TO (110,140) ,LAYNUM 

C#** KX 
Cc THIS SECTION IS FOR THE ONE LAYER CASE 
Cc 45874 ES VELOCI£Y IN LOWEST OF UPPER LAYERS 
C#* 

110 Vov=4.874/V 
CK eK 
¢ CALCULATE A FIRST APPROXIMATION TO ALPHA(1,1) 
C# KK 


FALPHA=X(1)/VT2 

DOM WShI=1532 

IF (FALPHA.LT.FUNAD(I1,J)) GO TO 116 
115 CONTINUE 
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STOP 11 
116 IALP=J-2-1ITHETA 
ALPHA (1,1) =IALP/57. 29578 
DELTA=0.0 
CK KX 
Cc LOOP ON I IS FOR NUMBER OF TRACES 
E LOOP ON J IS FOR 10 ITERATIONS IN AN ATTEMPT TO FIND 
CK KX RAY PATH 
DO 130 I=1,NUMTR 
DO 120 J=1,10 
S=VOV*SIN (ALPHA (1,1) ) 
ALPHA 1=ARSIN ($) 
BETA (1,1) =ALPHA (1,1) +2. 0*THETA 
$=VOV*SIN (BETA (1,1) ) 
BETA1=ARSIN($) 
ALPH1D=57.29578*ALPHA1 
ITALP1=ALPH1D 
I$=IALP1+31 
$=ALPH1D-IALP1 
INDEX=I$ 
IF (ALPHA1.LT.0.0) INDEX=INDEX-1 
IF (INDEX.GT.60.OR.INDEX.LT.1) STOP 12 
XP=UPFUN (14) +$* (UPFUN (INDEX+1) -UPFUN (INDEX) ) 
I2=IALP+1+ITHETA 
INDEXS=I12 
IF (ALPHA (1,1).LT.0.0) INDEXS=I2-1 
IF (INDEXS.GT.30.OR.INDEXS.LT.1) STOP 13 
AITDIV=(FUNAD(I1,INDEXS+1) -FUNAD (11, INDEXS) ) 
FUNAD1=FUNAD (11,12) +DELTA*AITDIV 
PATHNI=VT2+XP*STH 
TIMENI=PATHNI/V 
XE=XP+FUNAD1*PATHNI 
BETA1D=BETA1*57.29578 
IBETA1=BETA1D 
I$=IBETA1+31 
S=BETA1D-IEBETA1 
INDEX=I$ 
IF (BETA1.LT.0.0) INDEX=INDEX-1 
IF (INDEX.GT.60.OR.INDEX.LT.1) STOP 14 
XP=XP+UPFUN (IS) +$* (UPFUN (INDEX+1) -UPFUN (INDEX) ) 
$=XP-X (1) 
IF (ABS($).LT.0.025) GO TO 125 
C** eX 
Cc CALCULATE A NEXT APPROXIMATION TO ALPHA (1,1) 
C#* HX 
ALPHA (1,1) =ALPHA (1,1) -$/XP*FUNAD1/AITDIV/57. 29578 
ITALP=ALPHA (1,1) *57.29578 
120 DELTA=ALPHA (1,1) *57.29578-IALP 


STOP 15 
C** KX 
Cc CALCULATE TCOR(I)=TIME CORRECTION ON I'TH TRACE AND 
Cc ALSO CALCULATE A FIRST APPROXIMATION TO ALPHA(1,1) 
Cc FOR NEXT I VALUE 
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125 TCOR (I) =TIMENI* (FUNAED (11,12) +DELTA* (FUNAED (11, INDEXS+ 
- 1) -FUNAED (11, INDEXS)) ) 
TCOR (I) =TCOR (I) +UPFUNT (IBETA1+31) +UPFUNT (IALP 1431) 
IF (I.EQ.NUMTR) GO TO 200 
ALPHA (1,141) =ALPHA (1,1) +(X (I+1) -X (I) ) /X (I) *FUNAD1/AITD 
.1V/57. 29578 
IALP=ALPHA (1,141) *57. 29578 
DELTA=ALPHA(1,1+1) *57.29578-LALP 
130 CONTINUE 
STOP 16 
Cc 
CCK ek 
Cc STATEMENT 140 STARTS THE 2 LAYER CASE 
C** KX 
Cc 
140 VOV=4.874/V1 
C**** 4,874 IS THE VELOCITY IN LOWEST OF UPPER LAYERS 
V10V=V1/V 
CK KK 
Cc CALCULATE A FIRST APPROXIMATICN TO ALPHA (2,1) 
C¥* ** 
APFACT=VT2/(VT2+VT21) 
FALPHA=X (1) /VZ2*APFACT 
Be 4890451431 
IF (FALPHA.LT.FUNAD(I1,J)) GO TO 146 
145 CONTINUE 
STOP 21 
146 IALP=J-2-ITHETA+ITHET1 
ALPHA (2,1) =IALP/57.29578 
DELTA=0.0 
CK KK 
Cc LOOP ON I IS FOR NUMBER OF TRACES 
Cc LOOP ON J IS FOR 10 ITERATIONS IN AN ATTEMPT TO FIND 
G RAY PATH 
C%* * KK 
DO 160 I=1,NUMTR 
DO 150 J=1,10 
$=V1OV*SIN (ALPHA (2,1)) 
ALPHA (1,1) =ARSIN($) -THETA1 
$=VOV*SIN (ALPHA (1,1)) 
ALPHA1=ARSIN ($) 
BETA(2,1) =ALPHA(2,1) +2.0* (THETA-THETA1) 
$=V10V*SIN (BETA (2,1)) 
BETA (1,1) =ARSIN ($) +THETA1 
$=VOV*SIN (BETA(1,1) ) 
BETA1=ARSIN($) 
ALPH1D=ALPHA1*57. 29578 
IALP1=ALPH1D 
I$=IALP1+31 
$=ALPH1D-IALP1 
INDEX=I$ 
IF (ALPHA1.LT.0.0) INDEX=INDEX-1 
IF (INDEX.GT.60.OR.INDEX.LT.1) STOP 22 
XP=UPFUN (I$) +$* (UPFUN (INDEX+1) -UPFUN (INDEX) ) 
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D1=VT21+XP*STH1 
D1=D1/COS (ALPHA (1,1) +THETA1) 

CK KK 

Cc D1 IS LENGTH OF RAY PATH IN LAYER 1 

C**K KK 
XP=XP+D1*SIN (ALPHA (1,1) ) 
I 2=IALP+1+ITHETA-ITHET1 
INDEXS=12 
IF (ALPHA (2,1I).LT.0.0) INDEXS=I2-1 
IF (INDEXS.GT.30.OR.INDEXS.LT.1) STOP 23 
AITDIV=FUNAD (11, INDEXS+1) -FUNAD (11, INDEXS) 
FUNAD1=FUNAD (11,12) +DELTA*AITDIV 
PATHNI=VT2+XP*STHM1/CTH1 
TIMENI=PATHNI/V 
XP=XP+FUNAD1*PATHNI*CTH1 
$=VT21/CTH14+XP*STH1/CTH1 
D1=D1+$/COS(BETA(1,1)) 
XP=XP+$*TAN (BETA (1,1) ) 
BETA1D=57. 29578*BETA1 
IBETA1=BETA1D 
I$=IBETA1+31 
$=BETA1D-IBETA1 
INDEX=I$ 
IF (BETA1.LT.0.0) INDEX=INDEX-1 
IF (INDEX.GT.60.OR.INDEX.LT.1) STOP 24 
XP=XP+UPFUN (IS) +$* (UPFUN (INDEX+1) -UPFUN (INDEX) ) 
$=XP-X (I) 
IF (ABS ($).LT.0.025) GO TO 155 

C#* *X 

Cc CALCULATE A NEXT APPROXIMATION TO ALPHA (2,1) 

CK OK 
ALPHA (2,1) =ALPHA (2, 1) -$/XP*¥FUNAD1/AITDIV/57. 29578 
IALP=ALPHA (2,1) *57.29578 

150 DELTA=ALPHA (2,1) *57.29578-IALP 


STOP S25 
CX KK 
C CALCULATE TCOR(1I)=TIME CORRECTION ON I'TH TRACE AND 
Cc ALSO CALCULATE A FIRST APPROSIMATION TO ALPHA(2,i) FOR 
c NEXT I VALUE 
C*K KK 


155 TCOR (I) =TIMENI* (FUNAED (11,12) +DELTA* (FUNAED (11, 1NDEXS 
.+1) -FUNAED (11, INDEXS) )) 
TCOR(L) =TCOR (I) +D1/V1+UPFUNT (IBETA1+31) tUPFUNT (IALP1+3 
ony 
IF (I.EQ.NUMTR) GO TO 200 
ALPHA (2,1+1) =ALPHA (2, 1) + (X (I+1) -X (I) ) /X (I) *FUNAD1/AITD 
. IV*APFACT/57. 29578 
ITAL P=ALPHA (2,1+1) *57. 29578 
DELTA=ALPHA(2,1+1) *57.29578-IALP 
160 CCNTINUE 
STOP 26 
200 DO 210 I=1,NGATE 
SUM (I) =0.0 
210 SUMNOR (i) =0.0 
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CALCULATE ISPEC AS A MEASURE OF GOODNESS OF FIT OF 
REFLECTION ARRIVALS TO A PARTICULAR STEPOUT 
SPECTELED BY ¥), (LHETA, tAND .TO 


DO 250 I=1,NUMTR 
INDEX=178.5714* (TCOR (I) -TIMINC) -IP+2.5 
DO 240 J=1,NGATE 

IF (INDEX.LT.1.0OR.INDEX.GT.500) STOP 5 
SUM (J) =SUM (J) +A (INDEX, I) 

SUMNOR (J) =SUMNOR (J) +ABS (A (INDEX,1)) 
INDEX=INDEX+1 

CONTINUE 

SUMSQ=0.0 

SUMNSQ=0.0 

DO 260 I=1,NGATE 
SUMSQ=SUMSQ+SUM (I) *SUM(Z) 
SUMNSQ=SUMNS Q+SUMNOR (1) *SUMNOR (1) 
ISPEC (ICT3)=0 


IF (SUMNSQ.NE.0.0) ISPEC(ICT3) =SUMSQ/SUMNSQ* 1000. 


V=V+VELINC 


C****WRITE TIME CORRECTIONS ON DISK FOR GOOD CORRELATIONS. 


Cc 


300 


400 


Ol 
ic 
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IF (ISPEC (ICT3).GT.800) WRITE (3) ICT1,ICcT 
(I) ,I=1,NUMTR) 

CONTINUE 

WRITE (6,7) TO, (ISPEC (I) ,I=1,NVEL) 

TO=TO4+TIMINC 

CONTINUE 

ITHETA=ITHETA+ITHINC 

CONTINUE 

FORMAT (2034) 

FORMAT (315, 2F5.2,15) 

FORMAT (12F5. 3) 

FORMAT (2F5.3,15,2F5.2,415) 

FORMAT ('1',20A4," AT A DIP OF ',13) 

FORMAT (*-TIME VEL=',20(F4.2,2X)) 

FORMAT (' ',F5.2,5X,20(14,2X)) 

FORMAT ('1X DISTANCES IN KM. FOR ',20A4,'° 
', 12F10. 3) 


9 FORMAT (*-LAYER NUMBER = ',11) 
10 FORMAT('OV1 = ',F5.3,5X,'T1 = ',F4.2,5X,"THETA1 = ',13 


) 


S99 StL 


END 


2,1CT3, (TCOR 


FOLLOW: ',//' 


ive 
Oa 

% 
on 
= 
r 


io TI 1 sexucsod a) sein dee B G&A 


3 
gate age orrea 4 OF SEAVI 


eX @TAIOOTAD 3 
an& WO IT! a 
fh Gi thas 


ATUOS, Par o2c og 
-2*2i~-({ WIATI- (1) 8ODT) PRP TEUSVE KK aL 
bg ad th vas Of 

(000.70 .xX% 40 ree *RGGNI) <I 
(t . X0 > (i) AUR (b) nv 


((1 4S) €) 2G + eb) LOAM VES (by SOWME ; : 
r+ Su I=3EGM One 
UMEPEOS O2S 


; o: .o=G8a02 


0-0 20eR, 


TAGE, THI OG OB 
(5) (L) HUe4+Oesi2e=Ge sug 
ha) 3 (A) TUC+QCHRGES =O CHMVESY 
Q= (ESOL) Da Tart 
- - “as Z 
2) - Le she CARVE) i 
bon * ; 
< SRLIZV+V=9 
‘ A { £VIOAaRHOD: BATT =z 
7 Zz : a! : 3 Oie we Fe C27 ji} Dass} pi 5 
=< tan) = 
= . (HTHOU,T = ‘ (1) < 
aad oe 


 * 
1LLHTESATHESTISAT ANTI 
: SUWITHOD Oe 
¥ (HAGE) Tanaot P © 
(2E{R. 29S ERE. Lan ao® & 
(tC -CERS). FAR IOS 
Lm 4 ¥ 


£x, = TAPEBT KEy6 09! 4 £T*) 
ad 


i od A -.. 
*. Oe = 


(CHEE) ey Phe wane 


82 


SUBROUTINE RGET (LUNIT,1IN,A,N,ICODE,IP,IC, *) 


THIS SUBROUTINE READS THE DATA OFF TAPE WITH THE AID OF 
THE SYSTEM SUBROUTINE READ. THE DATA WAS WRITTEN ON 
THE TAPE AS VS TYPE RECORDS WITH BLOCKSIZE OF 2*IN-8. 
THE BLOCKSIZE OF THE TAPE MUST BE < OR = 8200. 

N IS THE NUMBER OF DATA POINTS DESIRED STARTING AT 
POINT IP IN IA. IC IS THE CHANNEL NUMBER OF THE 

CHANNEL THAT IS TO BE SKIPPED. IF IC=0 NO CHANNEL IS 
SKIPPED. ICODE IS THE NUMBER OF CHANNELS TO BE READ. 


INTEGER¥2 INLEN 

REAL IA (4096) 

DIMENSION A(1) 

INLEN=IN 

IEND=ICODE 

FORMAT ('-BAD RETURN FROM READ IN GET') 
K=0 

IF (IC.NE.O) IEND=IEND+1 

DO 3 I=1,1END 

CALL READ (IA,INLEN,0, LNR, LUNIT,&4) 
T Petre RO.i1C)) Go Tas 

DOrsed=17N 

A (K+J) =IA (IP +J-1) 

K=K+N 

CONTINUE 

RETURN 

WRITE (6,1) 

RETURN 1 

END 
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